/*
*  multipole.cpp - routines to perform multipole calculations.
*
*  W. Elliott
*  W. T. Rankin
*
*  Copyright (c) 1994 Duke University
*     All rights reserved.					*	*/

/*			notes on compress version..
			we use smallcell with addmultipoleC
			Calc_M2M parent.Smallcell -child.Smallcell

*/

#include "multipole.h"
#include "mpole_misc.h"



/****************************************************************
*
*  this procedure performs the one-time mpe initialization of
*  global constants and arrays
*
*/

int MultipoleSetup( int Mp )
{
	//cout<<"in multipole setup Mtermp = " << Mtermp <<endl;
   Cinit(Mp);
   return 0;
}

/*****************************************************************
*
*  AddMultipoleC calculates the Coulomb multipole expansion of a particle 
*  of charge q, at location <vect> relative to the origin. The result is
*  accumulated into the blank multipole expansion array M. The array size
*  is p.
*
*/

void AddMultipoleC( Kompleks **M, int p, Real q, vektor vect)
{
   int   n, m;
   vektor sv;
   /*
    * convert input to a spherical vector
    */
   Cart2Sph(vect, &sv);

   /*
    * compute the F matrix, leaving the result in Y_C{n][m]
    */

   makeF(p, sv);

   for (n = 0; n < p; n++) {
      for (m = 0; m <= n; m++) {
	 M[n][m].x += q * Y_C[n][m].x;
	 M[n][m].y -= q * Y_C[n][m].y;
      }	/* for m */
   } /* for n */
} /* MakeMultipoleC */

/****************************************************************
*
*  Calc_M2M() - shift the multipole expansion for a child cell and add
*    it to the multipole expansion of the parent cell.
*
*  w.t.rankin
*
*/
int Calc_M2M( Node* ccell, Node* pcell)
/* child cell */  /* parent cell */
   {
   Mtype cmpole;      /* child multipole expansion */
   Mtype pmpole;      /* parent multipole expansion */
   vektor v;          /* position vector from parent to child center */

   cmpole = ccell->m;
   pmpole = pcell->m;

   /*
   *  calculate the cartesian vector from the center of the child cell
   *  to the center of the parent cell
   */

   v = pcell->SmallCell - ccell->SmallCell;

   /*
   *  shift the child multipole expansion, placing the result in 
   *  the parent multipole expansion
   */

   //cout<<"-- M2M olan child ="<<ccell->level<<" parent="<<pcell->level;

   M2M_C(cmpole, pmpole, Mtermp, v);

   return 0;
} 

/****************************************************************
*
*  Calc_M2L() - shift the multipole expansion for a remote cell and add
*    it to the local expansion of the local cell.
*
*  w.t.rankin
*
*/

int Calc_M2L(Node *rcell, Node *lcell)
//int rcell;        /* remote cell */
//int lcell;        /* local cell */
   {

   Mtype rmpole;      /* remote multipole expansion */
   Mtype llocal;      /* local cells local expansion */
   vektor v;          /* position vector from parent to child center */


   /*
   *  check to see if the remote cell has a valid multipole expansion
   *   
   if ( rcell->mvalid == false )
   {
	   cout<<"M2L icinde false var."<<endl;
      return -1;
   }
   //else cout<<"M2L icinde true var."<<endl;

   /*
   *  set the local multipole expansion to valid
   */

   //lcell->lvalid = true;

   rmpole = rcell->m;
   llocal = lcell->l;

   /*
   *  calculate the cartesian vector from the center of the remote cell
   *  to the center of the local cell
   */

   //v = lcell->geocenter - rcell->geocenter;
   //v = lcell->geocenter - rcell->SmallCell;
   v = lcell->SmallCell - rcell->SmallCell;
   
   /*
   *  shift the remote multipole expansion, placing the result in 
   *  the local cells local expansion
   */

   M2L_C(rmpole, llocal, Mtermp, v);

   return 0;

   } 



/****************************************************************
*
*  Calc_L2L() - take the local expansion for a parent cell and
*     add it to the local expansion of a child cell.
*
*  w.t.rankin
*
*/

void Calc_L2L( Node* pcell, Node* ccell)
    /* parent cell */  /* child cell */
{

   Mtype plocal;      /* parent local expansion */
   Mtype clocal;      /* child local expansion */
   vektor v;          /* position vector from parent to child center */


   /*
   *  check to see if the parent cell has a valid local expansion
   * 
   if ( pcell->lvalid == false )
   {
	   //cout<<"L2L icinde lvalid false var parent "<<pcell->index <<" level = " << pcell->level<<endl;
	   //Display_local( pcell );
      return;
   }

   /*
   *  set the parent multipole expansion to valid
   */

   //ccell->lvalid = true;

   plocal = pcell->l;
   clocal = ccell->l;

   /*
   *  calculate the cartesian vector from the center of the parent cell
   *  to the center of the child cell
   */

   //v = ccell->geocenter - pcell->geocenter;
   v = ccell->SmallCell - pcell->SmallCell;
   
   /*
   *  shift the remote multipole expansion, placing the result in 
   *  the local multipole expansion
   */

   L2L_C(plocal, clocal, Mtermp, v);

} 



/****************************************************************
*
* misc routines to handle multipoles
*
*/

void Clear_mpole( Node* thiscell)
{
   /* 
   *  zero the expansion 
   */

   CMclear( thiscell->m , Mtermp);

   //thiscell->mvalid = false;
}

void Clear_local( Node* thiscell)
{
   /* 
   *  zero the expansion 
   */

   CMclear(thiscell->l ,Mtermp);

   //thiscell->lvalid = false;

}

/****************************************************************
*
*  CMclear(M,p) -
*     clears the multipole expansion of size p pointer to by Mptr
*
*/

int CMclear( Kompleks **M , int p )
{
   int i,j;   /* loop counters */

   for ( i=0; i<p; i++) {
      for ( j=0; j<=i; j++ ) {
         M[i][j].x = 0.0;
         M[i][j].y = 0.0; 
         }
      }
   return 0;
}



void display_expansion ( Kompleks **M , int p )
{
   int i,j;   /* loop counters */

   for ( i=0; i<p; i++) {
      for ( j=0; j<=i; j++ ) {
         cout<<"M["<<i<<"]["<<j<<"].x="<<M[i][j].x;
		 cout<<" - M["<<i<<"]["<<j<<"].y="<<M[i][j].y;
		 cout<<endl;
         }
      }
   cout<<"********************************************"<<endl;
}


void Display_mpole( Node* thiscell )
{
	cout<<"m of Node ... "<<thiscell->index<<endl;
	display_expansion( thiscell->m , Mtermp ); 
}

void Display_local( Node* thiscell )
{
	cout<<"l of Node ... "<<thiscell->index<<endl;
	display_expansion( thiscell->l , Mtermp ); 
}

/*****************************************************************
*
*  allocate and initialize static data structures needed by the
*  multipole calculations.
*
*/

void Cinit(int p )
{
   int      i, m, n;
   Real     *scratch, *factorial;
   Kompleks  *c_scratch;

   scratch = (Real *) malloc(((p * (p + 1)) / 2) * sizeof(Real));
   LegPoly = (Real **) malloc(p * sizeof(Real *));
   for (i = 0; i < p; i++) {
      LegPoly[i] = scratch;
      scratch += i + 1;
   } /* for i */

   scratch = &LegPoly[0][0];
   for (n = 0; n < (p * (p + 1)) / 2; n++) {
      scratch[n] = 0.0;
   }

   c_scratch = (Kompleks *) malloc(((p * (p + 1)) / 2) * sizeof(Kompleks));
   Y_C = (Mtype ) malloc(p * sizeof(Kompleks *));
   for (n = 0; n < p; n++) {
      Y_C[n] = c_scratch;
      c_scratch += n + 1;
   } /* for n */

   scratch = &Y_C[0][0].x;
   for (n = 0; n < (p * (p + 1)); n++)
      scratch[n] = 0.0;

   c_scratch = (Kompleks *) malloc(((p * (p + 1)) / 2) * sizeof(Kompleks));
   L = (Mtype ) malloc(p * sizeof(Kompleks *));
   for (n = 0; n < p; n++) {
      L[n] = c_scratch;
      c_scratch += n + 1;
   } /* for n */

   scratch = &L[0][0].x;
   for (n = 0; n < (p * (p + 1)); n++) {
      scratch[n] = 0.0;
   }

   Yxy = (Kompleks *) malloc((p + 1) * sizeof(Kompleks));

   factorial = (Real *) malloc(2 * (p + 1) * sizeof(Real));
   factorial[0] = 1.0;
   for (n = 1; n < 2 * (p + 1); n++) {
      factorial[n] = (Real) n *factorial[n - 1];
   } /* for n */

   scratch = (Real *) malloc(((p * (p + 1)) / 2) * sizeof(Real));
   Fcoeff = (Real **) malloc(p * sizeof(Real **));
   for (n = 0; n < p; n++) {
      Fcoeff[n] = scratch;
      scratch += n + 1;
   } /* for n */

   for (n = 0; n < p; n++) {
      for (m = 0; m <= n; m++) {
    Fcoeff[n][m] = pow(-1.0, (Real) (n + m)) / factorial[n + m];
      } /* for m */
   } /* for n */

   scratch = (Real *) malloc(((p * (p + 1)) / 2) * sizeof(Real));
   Gcoeff = (Real **) malloc(p * sizeof(Real **));
   for (n = 0; n < p; n++) {
      Gcoeff[n] = scratch;
      scratch += n + 1;
   } /* for n */

   for (n = 0; n < p; n++) {
      for (m = 0; m <= n; m++) {
    Gcoeff[n][m] = pow(-1.0, (Real) (n + m)) * factorial[n - m];
      } /* for m */
   } /* for n */

   free(factorial);

} /* Cinit */

/****************************************************************
 *
 * M2M_C translates a Coulomb multipole expansion M1 to another
 * expansion M2, both size p, along the vektor v.  the vektor points
 * from the center of M1 to the center of M2.
 *
 */
void M2M_C(Kompleks **M1, Kompleks **M2 , int p , vektor vect )
{
   int             n, m, np, mp, startm, endm;
   Real            atemp;
   vektor          sv;

   /*
    * convert input to a spherical vektor
    */

   Cart2Sph(vect,&sv);

   /*
    * compute the F matrix, leaving the result in Y_C{n][m]
    */

   sv[0] = -sv[0];
   makeF(p, sv);

   for (np = 0; np < p; np++) {
      for (mp = 0; mp <= np; mp++) {
    for (n = 0; n <= np; n++) {
       startm = mp - (np - n);
       startm = (startm < -n ? -n : startm);
       endm = mp + (np - n);
       endm = (endm > n ? n : endm);
       if (startm <= endm) {
          m = startm;
          while (m < 0 && m <= endm) {
        atemp = 1.0 - 2.0 * (Real) (0x0001 & (-m));
        /* above eqn = -1^m */
        M2[np][mp].x += atemp *
           (M1[n][-m].x * Y_C[np - n][mp - m].x -
            M1[n][-m].y * Y_C[np - n][mp - m].y);
        M2[np][mp].y -= atemp *
           (M1[n][-m].x * Y_C[np - n][mp - m].y +
            M1[n][-m].y * Y_C[np - n][mp - m].x);
        m++;
          }      /* while m negative */
          while (m < mp && m <= endm) {
        M2[np][mp].x +=
           (M1[n][m].x * Y_C[np - n][mp - m].x +
            M1[n][m].y * Y_C[np - n][mp - m].y);
        M2[np][mp].y +=
           (M1[n][m].x * -Y_C[np - n][mp - m].y +
            M1[n][m].y * Y_C[np - n][mp - m].x);
        m++;
          } /* while m, positive, less than mp */
          while (m <= endm) {
        atemp = 1.0 - 2.0 * (Real) (0x0001 & (mp + m));
        M2[np][mp].x += atemp *
           (M1[n][m].x * Y_C[np - n][m - mp].x -
            M1[n][m].y * Y_C[np - n][m - mp].y);
        M2[np][mp].y += atemp *
           (M1[n][m].x * Y_C[np - n][m - mp].y +
            M1[n][m].y * Y_C[np - n][m - mp].x);
        m++;
          } /* while m */
       } /* if startm <= endm (ie m exists) */
    } /* for n */
      } /* for mp */
   } /* for np */

} /* M2M_C */

/*****************************************************************
*
*  M2L_C converts a Coulomb multipole expansion M in to a local 
*  expansion L, shifting the original expansion along the vektor <vect>.
*  Both expansions are of size p.
*
*/

int M2L_C( Kompleks **M, Kompleks **L , int p , vektor vect )
{
   int n, m, np, mp;
   Real atemp,btemp;
   kompleks         *Mp, *Lp, *Hp;
   vektor sv;

   Cart2Sph(vect, &sv);
   /*
    * compute the G matrix, leaving the result in Y_C[n][m]
    */

   makeG(p, sv);


   Lp = &(L[0][0]);
   for (np = 0; np < p; np++) {
      atemp = 1.0;
      for (mp = 0; mp <= np; mp++) {
	 for (n = 0; n < p - np; n++) {
	    Mp = &(M[n][n]);
	    Hp = &(Y_C[np+n][n-mp]);
	    for (m = -n; m < -mp; m++) {
	       Lp->x += atemp * (Mp->x * Hp->x - Mp->y * Hp->y);
	       Lp->y -= atemp * (Mp->x * Hp->y + Mp->y * Hp->x);
	       Mp--;
	       Hp--;
	    } /* for m, m neg and |m| > |mp| */
	    Hp = &(Y_C[np+n][mp+m]);
            btemp = 1.0 - 2.0 * (Real) (0x0001 & (-m));
	    for (; m < 0; m++) {
	       Lp->x += btemp * (Mp->x * Hp->x + Mp->y * Hp->y);
	       Lp->y += btemp * (Mp->x * Hp->y - Mp->y * Hp->x);
	       Mp--;
	       Hp++;
	       btemp = -btemp;
            } /* for m, m neg and |m| <= |mp| */
	    for (; m <= n; m++) {
	       Lp->x += (Mp->x * Hp->x - Mp->y * Hp->y);
	       Lp->y += (Mp->x * Hp->y + Mp->y * Hp->x);
	       Mp++;
	       Hp++;
	    } /* for m, m neg and mp-m positive */
	 } /* for n */

	 Lp++;
	 atemp = -atemp;

      } /* for mp */
   } /* for np */

   return 0;

}  /* M2L_C */

/*****************************************************************
*
*  L2L_C translates a local expansion L1 to another local expansion L2, 
*  both of size p, along vektor <vect>
*
*/

int L2L_C( Kompleks **L1, Kompleks**L2 , int p , vektor vect ) 
{
   int np, mp, n, m;
   int startm, endm;
   Real atemp;
   vektor sv;

   Cart2Sph(vect, &sv);
   /*
    * compute the F matrix, leaving the result in Y_C{n][m]
    */

   sv[0] = -sv[0];
   makeF(p, sv);


   for (np = 0; np < p; np++) {
      for (mp = 0; mp <= np; mp++) {
	 for (n = np; n < p; n++) {
	    startm = mp - (n - np);
	    startm = (startm < -n ? -n : startm);
	    endm = mp + (n - np);
	    endm = (endm > n ? n : endm);
	    m = startm;
	    if (startm <= endm) {
	       while (m < 0 && m <= endm) {
		  atemp = 1.0 - 2.0 * (Real) (0x0001 & (mp));
		  L2[np][mp].x += atemp *
		     (L1[n][-m].x * Y_C[n - np][-m + mp].x +
		      L1[n][-m].y * Y_C[n - np][-m + mp].y);
		  L2[np][mp].y += atemp *
		     (L1[n][-m].x * Y_C[n - np][-m + mp].y -
		      L1[n][-m].y * Y_C[n - np][-m + mp].x);
		  m++;
	       } /* while m is negative */
	       while (m < mp && m <= endm) {
		  atemp = 1.0 - 2.0 * (Real) (0x0001 & (m + mp));
		  L2[np][mp].x += atemp *
		     (L1[n][m].x * Y_C[n - np][-m + mp].x -
		      L1[n][m].y * Y_C[n - np][-m + mp].y);
		  L2[np][mp].y += atemp *
		     (L1[n][m].x * Y_C[n - np][-m + mp].y +
		      L1[n][m].y * Y_C[n - np][-m + mp].x);
		  m++;
	       } /* while m-mp is negative */
	       while (m <= endm) {
		  L2[np][mp].x +=
		     (L1[n][m].x * Y_C[n - np][m - mp].x +
		      L1[n][m].y * Y_C[n - np][m - mp].y);
		  L2[np][mp].y +=
		     (L1[n][m].x * -Y_C[n - np][m - mp].y +
		      L1[n][m].y * Y_C[n - np][m - mp].x);
		  m++;
	       } /* while mp < endm, m and m-mp positive */
	    } /* if startm <= endm (ie m exists) */
	 } /* for n */
      } /* for mp */
   } /* for np */

   return 0; 

} /* L2L_C */


/*****************************************************************
*
*  Force_C computes the Coulomb force and potential on the particle with
*  charge q at position <vect>, due to the local expansion L of size p. 
*  The results are returned as a potential 'rnt_pot', and as a vektor
*  force 'rnt_force'.
*
*/

void Force_C(Kompleks **Lin, int p, Real q, vektor pv, double *potp, vektor *fv) 
{
   int             n, m;
   Real            pot, fr, fa, fb;
   Real            cosalpha, sinalpha, cosbeta, sinbeta;
   Real            cotalpha;
   Real            rtemp, ntemp, mtemp, ftemp, ptemp;
   vektor       sv;



   /*
    * we used to convert the input to a spherical vektor
    *
    * but now we can replace all of this stuff with divisions!
    */

   Real x, y, z, xy, xy2;
   Real rho, r2, alpha, beta;

   x = pv[0]; y = pv[1]; z = pv[2];

   xy2 = x*x + y*y;
   xy = sqrt(xy2);
   r2 = xy2 + z*z;
   rho = sqrt(r2);

   if (rho == 0.0) {
      alpha = 0.0;
      cosalpha = 1.0;
      sinalpha = 0.0;
   }
   else {
      sinalpha = xy / rho;
      if (sinalpha < SMALL_THETA) {
	 sinalpha = 0.0;
	 if (z >= 0.0) {
	    cosalpha = 1.0;
	    alpha = 0.0;
	 }
	 else {
	    cosalpha = -1.0;
	    alpha = M_PI_2;
	 }

      }	
      else {
	 cosalpha = z / rho;
	 alpha = acos(cosalpha);
      }
   }
   cotalpha = cosalpha/sinalpha;

   if ((x == 0.0) && (y == 0.0)) {
      cosbeta = 1.0;
      sinbeta = 0.0;
      beta = 0.0;
   }
   else {
      cosbeta = x / xy;
      sinbeta = y / xy;
      beta = atan2(y, x);
   }

   sv[0] = rho; sv[1] = alpha; sv[2] = beta;


   pot = 0.0;
   fr = 0.0;
   fa = 0.0;
   fb = 0.0;

   /*
    * compute the F matrix, leaving the result in Y_C{n][m]
    */

   sv[0] = -sv[0];
   makeF(p, sv);


   /*
    * find the potential by computing the product of L[n,m].F*[n,m]
    */

   for (n = 0; n < p; n++) {
      for (m = 0; m<=n; m++) {
         if (m == 0) {
            pot += Y_C[n][m].x*Lin[n][m].x + Y_C[n][m].y*Lin[n][m].y;
	 } /* if */
         else {
            pot += 2.0 * (Y_C[n][m].x*Lin[n][m].x + Y_C[n][m].y*Lin[n][m].y);
	 } /* else */
      } /* for m */
   }/* for n */


   /*
    *  this section computes the special case for the
    *  condition where either alpha=0 (assumes beta=0)
    *  or rho=0 (assume alpha=0 and beta=0)
    *
    */

   if ( sinalpha == 0.0 ) {

      if (sv[0] == 0.0) {

         fr += Fcoeff[1][0] * LegPoly[1][0] * Lin[1][0].x;
         fr += 2.0 * Fcoeff[1][1] * LegPoly[1][1] * Lin[1][1].x;

         /*
	  * rotate the r-vektor to the x-axis and recompute
	  * we can probably get away with only computing the
	  * first couple terms in the Legerdre polynomial
	  */

	 Legendre(LegPoly, p, 0.0);

         fa += Fcoeff[1][0] * LegPoly[1][0] * Lin[1][0].x;
         fa += 2.0 * Fcoeff[1][1] * LegPoly[1][1] * Lin[1][1].x;

         fb += Fcoeff[1][0] * LegPoly[1][0] * Lin[1][0].y;
         fb += 2.0 * Fcoeff[1][1] * LegPoly[1][1] * Lin[1][1].y;

      } /* sv[0] == 0 */


      /*
       * otherwise, we can compute the vektors in a more simple
       * fashion.
       */

      else {
         rtemp = 1.0 / sv[0];
	 ftemp = rtemp;
         for (n = 1; n < p; n++) {
            fr += ftemp * (Y_C[n][0].x*Lin[n][0].x + Y_C[n][0].y*Lin[n][0].y);
            for (m = 1; m<=n; m++) {
               fr += 2.0 * ftemp *
		  (Y_C[n][m].x*Lin[n][m].x + Y_C[n][m].y*Lin[n][m].y);
            } /* for m */
	    ftemp += rtemp;
         } /* for n */

         /*
          *  compute the alpha and beta component of the gradient
          */

         rtemp = 1.0;
         ntemp = 0.0;
         for (n = 1; n < p; n++) {
            ntemp -= (double)n;
            ftemp  = 2.0 * Fcoeff[n][1] * rtemp * ntemp;
            fa += ftemp * Lin[n][1].x;
            fb += ftemp * Lin[n][1].y;
            rtemp *= sv[0];
         } /* for n */

      } /* else */

      *potp = q * pot;
      (*fv)[0] = q * fa;
      (*fv)[1] = q * fb;
      (*fv)[2] = q * fr;

   } /* alpha == 0 */


   /*
    * if alpha is not equal to zero, we can compute the vektors
    * in the traditional manner.
    *
    */

   else {

      rtemp = 1.0 / sv[0];
      ftemp = rtemp;
      for (n = 1; n < p; n++) {
         fr += ftemp * (Y_C[n][0].x*Lin[n][0].x + Y_C[n][0].y*Lin[n][0].y);
         for (m = 1; m<=n; m++) {
            fr += (ftemp+ftemp) *
	       (Y_C[n][m].x*Lin[n][m].x + Y_C[n][m].y*Lin[n][m].y);
         } /* for m */
         ftemp += rtemp;
      } /* for n */


      for (n = 0; n < p; n++) {
         ftemp = 2.0;
         for (m = 1; m<=n; m++) {
            fb -= ftemp * (Y_C[n][m].y*Lin[n][m].x - Y_C[n][m].x*Lin[n][m].y);
            ftemp += 2.0;
         } /* for m */
      }/* for n */
      fb /= (sv[0] * sinalpha);


      /*
       *  compute the alpha component of the gradient
       *  this is an ugly loop with lots of non-obvious
       *  temp vars, but it does use the least number of
       *  FP multiplies.
       */

      ntemp = 2.0;
      for (n = 1; n < p; n++) {

         fa -= ntemp *
	    ((Y_C[n][1].x*cosbeta + Y_C[n][1].y*sinbeta) * Lin[n][0].x +
             (Y_C[n][1].y*cosbeta - Y_C[n][1].x*sinbeta) * Lin[n][0].y );

         mtemp = 1.0;
         for (m = 1; m < n; m++) {

	    ftemp = 2.0 * mtemp * cotalpha;
	    ptemp = 2.0 * (ntemp+mtemp);

            fa +=  ftemp *
               (Y_C[n][m].x * Lin[n][m].x + Y_C[n][m].y * Lin[n][m].y);

            fa -= ptemp *
	       ((Y_C[n][m+1].x*cosbeta + Y_C[n][m+1].y*sinbeta) * Lin[n][m].x +
                (Y_C[n][m+1].y*cosbeta - Y_C[n][m+1].x*sinbeta) * Lin[n][m].y );

            mtemp += 1.0;
         } /* for m */

         fa += 2.0 * mtemp * cotalpha * 
            (Y_C[n][n].x*Lin[n][n].x + Y_C[n][n].y*Lin[n][n].y); 

         ntemp += 1.0;

      } /* for n */

      fa /= sv[0];

      /*
       *  rotate the spherical vektor coordinates to cartesian
       */

      *potp = q * pot;
      (*fv)[0]= q * (fr*sinalpha*cosbeta + fa*cosalpha*cosbeta - fb*sinbeta);
      (*fv)[1] = q * (fr*sinalpha*sinbeta + fa*cosalpha*sinbeta + fb*cosbeta);
      (*fv)[2] = q * (fr*cosalpha - fa*sinalpha);

   } /* else */

} /* Force_C */

/****************************************************************
 *
 *  ForceM_C() -
 *
 *  computes the Coulomb force and potential on the particle
 *  with charge q at position pv.[x,y,z], due to the Multipole
 *  expansion Min of size p.  The results are returned in potp and
 *  fv.[x,y,z]
 *
 */

void ForceM_C( Kompleks  **Min, int p, Real q, vektor  pv, Real *potp, vektor  *fv )
{
   int             n, m;
   Real            pot, fr, fa, fb;
   Real            cosalpha, sinalpha, cosbeta, sinbeta;
   Real            cotalpha;
   Real            rtemp, ntemp, mtemp, gtemp, ptemp;
   vektor       sv;

   /*
    * convert input to a spherical vektor
    *
    * we could optimize this section the same as Force_C but
    * DPMTA never uses this procedure.  maybe someday.
    *
    */

   Cart2Sph(pv,&sv);

   cosalpha = cos(sv[1]);
   sinalpha = sin(sv[1]);
   if (sinalpha < SMALL_THETA) {
      sinalpha = 0.0;
      if (cosalpha > 0.0)
    cosalpha = 1.0;
      else
    cosalpha = -1.0;
   }            /* if a small alpha */
   cosbeta = cos(sv[2]);
   sinbeta = sin(sv[2]);
   cotalpha = cosalpha/sinalpha;

   pot = 0.0;
   fr = 0.0;
   fa = 0.0;
   fb = 0.0;

   /*
    * compute the G matrix, leaving the result in Y_C{n][m]
    */

   makeG(p, sv);


   /*
    * find the potential by computing the product of M[n,m].G[n,m]
    */

   for (n = 0; n < p; n++) {
      pot += Y_C[n][0].x*Min[n][0].x - Y_C[n][0].y*Min[n][0].y;
      for (m = 1; m<=n; m++) {
         pot += 2.0 * (Y_C[n][m].x*Min[n][m].x - Y_C[n][m].y*Min[n][m].y);
      } /* for m */
   }/* for n */


   /*
    *  this section computes the special case for the
    *  condition where alpha=0 (assume beta=0)
    *  rho cannot equal zero because the multipole expansion does 
    *  not converge.
    */

   if ( sinalpha == 0.0 ) {

      rtemp = 1.0 / sv[0];
      gtemp = rtemp;
      for (n = 0; n < p; n++) {
         fr += gtemp * (Y_C[n][0].x*Min[n][0].x - Y_C[n][0].y*Min[n][0].y);
         for (m = 1; m<=n; m++) {
            fr += 2.0 * gtemp *
          (Y_C[n][m].x*Min[n][m].x - Y_C[n][m].y*Min[n][m].y);
         } /* for m */
         gtemp += rtemp;
      } /* for n */

      /*
       *  compute the alpha and beta component of the gradient
       */

      rtemp = 1.0 / sv[0];
      ntemp = 0.0;
      for (n = 1; n < p; n++) {
         rtemp /= sv[0];
         ntemp += (double)n;
         gtemp  = 2.0 * Gcoeff[n][1] * rtemp * ntemp;
         fa += gtemp * Min[n][1].x;
         fb -= gtemp * Min[n][1].y;
      } /* for n */


      fa /= sv[0];
      fb /= sv[0];

      *potp = q * pot;
      (*fv)[0] = q * fa;
      (*fv)[1] = q * fb;
      (*fv)[2] = q * fr;

   } /* alpha == 0 */


   /*
    * if alpha is not equal to zero, we can compute the vektors
    * in the traditional manner.
    *
    */

   else {

      /* compute the rho component of the gradient */

      rtemp = 1.0 / sv[0];
      gtemp = rtemp;
      for (n = 0; n < p; n++) {
         fr += gtemp * (Y_C[n][0].x*Min[n][0].x - Y_C[n][0].y*Min[n][0].y);
         for (m = 1; m<=n; m++) {
            fr += 2.0 * gtemp *
          (Y_C[n][m].x*Min[n][m].x - Y_C[n][m].y*Min[n][m].y);
         } /* for m */
         gtemp += rtemp;
      } /* for n */


      /* compute the beta component of the gradient */

      for (n = 1; n < p; n++) {
         gtemp = 2.0;
         for (m = 1; m<=n; m++) {
            fb += gtemp * (Y_C[n][m].y*Min[n][m].x + Y_C[n][m].x*Min[n][m].y);
            gtemp += 2.0;
         } /* for m */
      }/* for n */
      fb /= (sv[0] * sinalpha);


      /* compute the alpha component of the gradient */

      ntemp = 1.0;
      for (n = 1; n < p; n++) {
         fa -= ntemp *
       ((Y_C[n][1].x*cosbeta + Y_C[n][1].y*sinbeta) * Min[n][0].x -
             (Y_C[n][1].y*cosbeta - Y_C[n][1].x*sinbeta) * Min[n][0].y );

         mtemp = 1.0;
         for (m = 1; m < n; m++) {

       gtemp = 2.0 * mtemp * cotalpha;
       ptemp = 2.0 * (ntemp-mtemp);

            fa += gtemp * (Y_C[n][m].x*Min[n][m].x - Y_C[n][m].y*Min[n][m].y); 

            fa -= ptemp *
          ((Y_C[n][m+1].x*cosbeta + Y_C[n][m+1].y*sinbeta) * Min[n][m].x -
                (Y_C[n][m+1].y*cosbeta - Y_C[n][m+1].x*sinbeta) * Min[n][m].y );

            mtemp += 1.0;
         } /* for m */

         fa += 2.0 * mtemp * cotalpha * 
            (Y_C[n][n].x*Min[n][n].x - Y_C[n][n].y*Min[n][n].y); 

      ntemp += 1.0;
      } /* for n */

      fa /= sv[0];
      fa *= -1.0;

      /*
       *  rotate the spherical vektor coordinates to cartesian
       */

      *potp = q * pot;
      (*fv)[0] = q * (fr*sinalpha*cosbeta + fa*cosalpha*cosbeta - fb*sinbeta);
      (*fv)[1] = q * (fr*sinalpha*sinbeta + fa*cosalpha*sinbeta + fb*cosbeta);
      (*fv)[2] = q * (fr*cosalpha - fa*sinalpha);

   } /* else */

} /* ForceM_C */


/****************************************************************
 *
 *  makeG() -
 *
 *  compute the G transfer matrix, using the equation
 *  from Elliott's dissertation.  the results are stored in the
 *  global array Y_C[n][m].
 *
 */
void makeG( int p, vektor sv )
{
   int             n, m;
   Real            rinv, rtemp, ytemp;

   Legendre(LegPoly, p, cos(sv[1]));
   Fourier_C(p, sv[2]);

   rinv = 1.0 / sv[0];
   rtemp = rinv;
   for (n = 0; n < p; n++) {
      for (m = 0; m <= n; m++) {
	 ytemp = rtemp * Gcoeff[n][m] * LegPoly[n][m];
	 Y_C[n][m].x = ytemp * Yxy[m].x;
	 Y_C[n][m].y = ytemp * Yxy[m].y;
      } /* for m */
      rtemp *= rinv;
   } /* for n */

} /* makeG */


/****************************************************************
 *
 *  makeF() - 
 *
 *  compute the F transfer matrix, as specified in 
 *  Elliott's dissertation.  the results are stored in the global
 *  array Y_C[n][m].
 *
 */
void makeF( int p, vektor sv )
{
   int     m, n;
   Real    rtemp, ytemp;

   Legendre(LegPoly, p, cos(sv[1]));
   Fourier_C(p, sv[2]);

   rtemp = 1.0;

   for (n = 0; n < p; n++) {
      for (m = 0; m <= n; m++) {
	 ytemp = rtemp * Fcoeff[n][m] * LegPoly[n][m];
	 Y_C[n][m].x = ytemp * Yxy[m].x;
	 Y_C[n][m].y = ytemp * Yxy[m].y;
      } /* for m */
      rtemp *= sv[0];//sv[0]
   } /* for n */
} /* makeF */


/****************************************************************
 *
 *  Fourier_C() -
 *
 *  compute the Fourier series components of angle b out to p+1 (0->p)
 *  terms.  results are left in the global array Yxy.
 *
 */
static void Fourier_C( int p, Real b)
{
   int m;

   if (Yxy == NULL) {
      fprintf(stderr, "Fourier called with null pointer to array\n");
      exit(0);
   }
   for (m=0; m <= p; m++) {
      Yxy[m].x = cos((Real) m * b);
      Yxy[m].y = sin((Real) m * b);
   }
} /* fourier */


/*****************************************************************
*
*  allocates the multipole expansion of size p leaving the resulting
*  matrix in Mptr.
*
*/
int Calloc(Kompleks ***Mptr , int p ) 
{
   int n;
   Real *scratch;
   Kompleks *scratchC;
   Kompleks **M;
  

   scratchC = (Kompleks *) malloc(((p * (p+1))/2) * sizeof(Kompleks));
   M = (Kompleks **) malloc( p * sizeof(Kompleks **));
   for (n = 0; n < p; n++) {
      M[n] = scratchC;
      scratchC += n+1;
   } /* for n */
   scratch = &M[0][0].x;
   for (n=0; n < (p * (p+1)); n++)
      scratch[n] = 0.0;
   *Mptr = M;

   return 0;
} /* Calloc */


/****************************************************************
 *
 *  copyG() - 
 *
 *  calculate the the G matrix and store the results in
 *  Yout.
 *
 */

void copyG( Kompleks **Yout, int p, vektor v )
{
   int			n, m;
   vektor       sv;


   /* convert input to a spherical vector */

   Cart2Sph(v,&sv);

   /*
    * compute the G matrix, leaving the result in Y_C{n][m]
    *
    * it would probably be better to use a memcopy() here
    *
    */


   makeG(p, sv);

   for (n = 0; n < p; n++) {
      for (m = 0; m <= n; m++) {
    Yout[n][m].x = Y_C[n][m].x;
    Yout[n][m].y = Y_C[n][m].y;
      }   /* for m */
   } /* for n */

} /* copyG */

/****************************************************************
 *
 *  copyF() - calculate the the F matrix and store the results in
 *  Yout.
 *
 */

void copyF( Kompleks **Yout, int p, vektor v )
{
   int             n, m;
   vektor       sv;


   /* convert input to a spherical vector */

   Cart2Sph(v,&sv);

   /*
    * compute the F matrix, leaving the result in Y_C{n][m]
    *
    * it would probably be better to use a memcopy() here
    *
    */

   sv[0] = -sv[0];
   makeF(p, sv);

   for (n = 0; n < p; n++) {
      for (m = 0; m <= n; m++) {
    Yout[n][m].x = Y_C[n][m].x;
    Yout[n][m].y = Y_C[n][m].y;
      }   /* for m */
   } /* for n */

} /* copyF */
